www.gusucode.com > 循环自相关函数工具箱源码程序 > matlab代做 修改 程序循环自相关函数工具箱/cyclostationary_toolbox/cyclic_3rd_order_cumulant_fast.m
function C3=cyclic_3rd_order_cumulant_fast(x1,x2,x3,T,max_tau) % % CYCLIC_3RD_ORDER_CUMULANT_FAST % % calculates the cyclic third order cumulant of % three signals x1,x2,x3 at frequency alpha using a fast % approximation based on the synchronous average of the % time varying autocorrelation. Fundamental signal % period can be defined as a single period or as a sequence % of once per period pulse times. % % C3(k*alpha,tau1,tau2)=E{(x1(t)-E{x1(t)}) * % (x2(t+tau1)-E{x2(t+tau1)} * % (x3(t+tau2)-E{x3(t+tau2)} * % exp(-jk(alpha)t) } % for k=0 ... 1/alpha % % % USAGE % C3=cyclic_3rd_order_cumulant_fast(x1,x2,x3,alpha,max_tau) % % File: cyclic_3rd_order_cumulant_fast.m % Last Revised: 25/11/97 % Created: 25/11/97 % Author: Andrew C. McCormick % (C) University of Strathclyde % Simple error checks if nargin~=5 error('Incorrect number of arguments for function cyclic_third_order_cumulant_fast'); end if T(1)<1 error('Synchronous period must be larger than 1 in function cyclic_third_order_cumulant_fast'); end [rows,cols]=size(x1); if rows>cols x1=x1'; end [rows,cols]=size(x2); if rows>cols x2=x2'; end [rows,cols]=size(x3); if rows>cols x3=x3'; end % Calculate synchronous averages from signals mx1=synchronous_average(x1,T); mx2=synchronous_average(x2,T); mx3=synchronous_average(x3,T); % Remove excess samples due to non-integer sampling % and renove cyclic mean from signal if length(T)==1 cp=1;np=1; while cp+T<length(x1) x1(cp:cp+floor(T)-1)=x1(cp:cp+floor(T)-1)-mx1; x2(cp:cp+floor(T)-1)=x2(cp:cp+floor(T)-1)-mx2; x3(cp:cp+floor(T)-1)=x3(cp:cp+floor(T)-1)-mx3; cp=cp+floor(T); np=np+T; if (np-cp)>1 x1=[x1(1:cp-1);x1(cp+1:length(x1))]; x2=[x2(1:cp-1);x2(cp+1:length(x2))]; x3=[x3(1:cp-1);x3(cp+1:length(x3))]; np=np-1; end end n=floor((length(x1)-2*max_tau-1)/T); else % extract time series correlated with periodic pulses rot_positions=T; T=floor(median(diff(rot_positions))); nx1=[]; nx2=[]; nx3=[]; n=length(rot_positions)-2; for k=1:n; cp=rot_positions(k); nx1=[nx1; x1(cp:cp+T-1)-mx1]; nx2=[nx2; x2(cp:cp+T-1)-mx2]; nx3=[nx3; x3(cp:cp+T-1)-mx3]; end nx1=[nx1;x1(rot_positions(n+1):rot_positions(n+1)+tau+1)-mx1(1:tau+1)]; x1=nx1; nx2=[nx2;x2(rot_positions(n+1):rot_positions(n+1)+tau+1)-mx2(1:tau+1)]; x2=nx2; nx3=[nx3;x3(rot_positions(n+1):rot_positions(n+1)+tau+1)-mx3(1:tau+1)]; x3=nx3; end % Compute time varying third order cumulant r=zeros(max_tau+1,max_tau+1,floor(T)); temp=zeros(floor(T),n); t=(1:floor(T)*n); for tau1=0:max_tau for tau2=0:max_tau temp(:)=x1(t).*x2(t+tau1).*x3(t+tau2); r(tau1+1,tau2+1,:)=mean(temp'); end end % Take DFT of time varying toc C3=zeros(max_tau+1,max_tau+1,floor(T)); for tau1=0:max_tau for tau2=0:max_tau C3(tau1+1,tau2+1,:)=fft(r(tau1+1,tau2+1,:))/T; end end